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£K Abstract 

| Fractal Interpolation has been proposed in the literature as an efficient way to construct clo- 

sure models for the numerical solution of coarse-grained Navier-Stokes equations. It is based on 
synthetically generating a scale-invariant subgrid-scale field and analytically evaluating its effects 
^ ' on large resolved scales. In this paper, we propose an extension of previous work by developing 

a multiaffine fractal interpolation scheme and demonstrate that it preserves not only the fractal 

O 

• <-h . dimension but also the higher-order structure functions and the non-Gaussian probability density 

function of the velocity increments. Extensive a-priori analyses of atmospheric boundary layer 

Oh' 

measurements further reveal that this Multiaffine closure model has the potential for satisfactory 
performance in large-eddy simulations. The pertinence of this newly proposed methodology in the 
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case of passive scalars is also discussed. 
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I. INTRODUCTION 



Generation of turbulence-like fields (also known as Synthetic Turbulence) has received 
considerable attention in recent years. Several schemes have been proposed 
with different degrees of success in reproducing various characteristics of turbulence. Re- 
cently, Scotti and Meneveau 0, [?| further broadened the scope of synthetic turbulence 
research by demonstrating its potential in computational modeling. Their innovative turbu- 
lence emulation scheme based on the Fractal Interpolation Technique (FIT) j^, |^] was found 
to be particularly amenable for a specific type of turbulence modeling, known as Large- 
Eddy Simulation (LES, at present the most efficient technique available for high Reynolds 
number flow simulations, in which the larger scales of motion are resolved explicitly and the 
smaller ones are modeled). The underlying idea was to explicitly reconstruct the subgrid 
(unresolved) scales from given resolved scale values (assuming computation grid-size falls in 
the inertial range of turbulence) using FIT and subsequently estimate the relevant subgrid- 
scale (SGS) tensors necessary for LES. Simplicity, straightforward extensibility for multi- 
dimensional cases, and low computational complexity (appropriate use of Fractal Calculus 
can even eliminate the computationally expensive explicit reconstruction step, see Section 
|lV|for details) makes this FIT-based approach an attractive candidate for SGS modeling in 
LES. 

Although the approach of ly, Ml is better suited for LES than any other similar scheme 

niJnnn 

(e.g., |2j, La, |4] , Ll ) , it falls short in preserving the essential small-scale properties of tur- 
bulence, such as multiaffinity (will be defined shortly) and non-Gaussian characteristics 
of the probability density function (pdf) of velocity increments. It is the purpose of this 
work to extend the approach of 0, Q| in terms of realistic turbulence-like signal generation 
with all the aforementioned desirable characteristics and demonstrate its potential for LES 
through a-priori analysis (an LES-SGS model evaluation framework). We will also demon- 
strate the competence of our scheme in the emulation of passive-scalar fields for which the 
non-Gaussian pdf and multiaffinity are significantly pronounced and cannot be ignored. 
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II. BASICS OF FRACTAL INTERPOLATION 



The fractal interpolation technique is an iterative affine mapping procedure to construct 
a synthetic deterministic small-scale field (in general fractal provided certain conditions 
are met, see below) given a few large-scale interpolating points (anchor points). For an 
excellent treatise on this subject, the reader is referred to the book by Barnsley |9]. In 
this paper, we will limit our discussion (without loss of generality) only to the case of three 
interpolating data points: {(xi,Ui) , i = 0,1,2}. For this case, the fractal interpolation 
iterative function system (IFS) is of the form {R 2 ; w n , n — 1, 2}, where, w n have the following 
affine transformation structure: 




a n 




n 



1,2. (1) 



To ensure continuity, the transformations are constrained by the given data points as fol- 

(xq \ I X n -\ \ I X 2 \ I x n \ 

= 1 and w n I = , for n = 1, 2. The parameters 

a n ,c n , e n and /„ can be easily determined in terms of d n (known as the vertical stretch- 
ing factors) and the given anchor points (xj,ttj) by solving a linear system of equations. 
The attractor of the above IFS, G, is the graph of a continuous function u : [xq, £2] — > R, 
which interpolates the data points (xj,Mj), provided the vertical stretching factors d n obey 
< \d n \ < 1. In other words, 

G = {(x, u (x)) : x E [x , x 2 ]} , 

where, (2) 
u (xi) = Ui, i = 0, 1, 2. 

Moreover, if \d\\ + |rf 2 | > 1 and ( not collinear, then the fractal (box-counting) 

dimension of G is the unique real solution D of |c?x| of -1 + \d- 2 \ a^~ l = 1 (for rigorous proof 
see 0). In the special case of three equally spaced points covering the unit interval [0,1], i.e., 
Xq = 0, x\ — 0.5 and x 2 = 1, the parameters of the affine transformation kernel become: 
a n = 0.5; c n = (u n - u n -i) - d n (u 2 - Uq) ; e n = x n _i, f n = u n _ x - d n u ; n = 1,2. In this 
case, the solution for the fractal dimension (D) becomes: 

D = l + log a (|d 1 | + |d 2 |) (3) 
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Notice that the scalings d% and d 2 are free parameters and cannot be determined using only 



, BE 



2j chose to use the 



equation (jSJ); at least one more constraint is necessary. For example, 
additional condition: \di\ = \d 2 \. 

III. SYNTHETIC TURBULENCE GENERATION 

Not long ago, it was found that turbulent velocity signals at high Reynolds numbers 
have a fractal dimension of D ~ 1.7 ± 0.05, very close to the value of D = 5/3 expected 
for Gaussian processes with a —5/3 spectral slope [lOj. For D = 5/3, the assumption 
of \di\ = \d 2 \ along with Equation El yields \di\ = \d 2 \ = 2 -1 / 3 Bill- One contribution 
of this paper is a robust way of estimating the stretching parameters without any ad-hoc 
prescription; the resulting synthetic field will not only preserve the fractal dimension (D) 
but also other fundamental properties of real turbulence. 

As an exploratory example, using the fractal interpolation IFS (Equation Q), we con- 
struct a 2 17 points long synthetic fractal series, u(x), with given coarse-grained points 



flQ: di 



(0.0, 1.2) , (0.5, —0.3) and (1.0, 0.7) and the stretching parameters used in 
-2- 1 / 3 ,d 2 = 2- 1 / 3 . Clearly, Fi gure la depicts that the synthetic series has fluctuations 
at all scales and it passes through all three interpolating points. 

Next, from this synthetic series we compute higher-order structure functions (see Figure 
lb for orders 2, 4 and 6), where the g^-order structure function, S q (r), is defined as follows: 

S q (r) = (\u(x + r) -u(x)\ q ) ~ r c « (4) 

where, the angular bracket denotes spatial averaging and r is a separation distance that varies 
in an appropriate scaling region (known as the inertial range in turbulence), 
exponent ( q is a nonlinear function of q, then following the convention of 
the field is called multiaffine, otherwise it is termed as monoaffine. In this context, we 
would like to mention that, Kolmogorov's celebrated 1941 hypothesis (a.k.a K41) based on 
the assumption of global scale invariance in the inertial range predicts that the structure 



If the scaling 

0, fl 3, fl I 



functions of order q scale with an exponent q/3 over inertial range separations U, l^ . 
Deviations from ( q = q/3 would suggest inertial range intermittency and invalidate the K41 



hypothesis. Inertial range intermittency is stil 



evidence for its existence is overwhelming [ill . 13^ . To interpret the curvilinear behavior of 



an unresolved issue, although experimental 



the ( q function observed in experimental measurements (e.g., Parisi and Frisch 
proposed the multifractal model, by replacing the global scale invariance with the assumption 
of local scale invariance. They conjectured that at very high Reynolds number, turbulent 
flows have singularities (almost) everywhere and showed that the singularity spectrum is 
related to the structure function-based scaling exponents, ( q by the Legendre transformation. 

Our numerical experiment with the stretching parameters of fill], i.e., \d\\ = jc^l — 2 -1 / 3 , 
revealed that the scaling exponents follow the K41 predictions (after ensemble averaging 
over one hundred realizations corresponding to different initial interpolating points), i.e., 
( q = q/3 (not shown here), a signature of monoaffine fields. Later on, we will give analytical 
proof that indeed this is the case for \di\ = |d 2 | = 2 -1 / 3 . Also, in this case, the pdfs of 
the velocity increments, 5u r (x) = u (x + r) — u (x), always portray near-Gaussian (slightly 
platykurtic) behavior irrespective of r (see Figure lc). This is contrary to the observations 
[111 ll 3l | . where, typically the pdfs of increments are found to be r dependent and become 
more and more non-Gaussian as r decreases. Theoretically, non-Gaussian characteristics of 
pdfs correspond to the presence of intermittency in the velocity increments and gradients 
(hence in the energy dissipation) 



121. 



In Figure 2, we plot the wavelet spectrum of this synthetic series. Due to the dyadic nature 
of the fractal interpolation technique, the Fourier spectrum will exhibit periodic modulation 

n 

(see Figures 7 and 8 of [7]). To circumvent this issue we make use of the (dyadic) discrete 



Haar wavelet transform. Following |15| . the wavelet power spectral density function E(K m ) 
is defined as: 

l{WT^ [i}) 2 )dx 



E (K m ) = — ~,'\ ' (5) 

v ; 27rln(2) v ' 

where, wavenumber K m {= i^z) corresponds to scale R m {= 2 m dx). The scale index m 

runs from 1 (finest scale) to log2(N) (coarsest scale). WT^ [i], dx and denote the Haar 

wavelet coefficient at scale m and location i, spacing in physical space and length of the 

spatial series, respectively. The power spectrum displays the inertial range slope of —5/3, 

as anticipated. 

At this point, we would like to invoke an interesting mathematical result regarding the 
scaling exponent spectrum, ( q , of the fractal interpolation IFS [ffil ]: 

JV 

C 9 = l-lo gjV ^K| 9 (6) 

n=l 
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where, N = the number of anchor points —1 (in our case N = 2). The original formulation 
was in terms of a more general scaling exponent spectrum, r (g) , rather than the 
structure function based spectrum ( q . The r (q) spectrum is an exact Legendre tranform 
of the singularity spectrum in the sense that it is valid for any order of moments (includ- 
ing negative) and any singularities 17, Q|. r (q) can be reliably estimated from data by 
the Wavelet-Transform Modulus-Maxima method \v\ . Q] . To derive Equation El from the 
original formulation, we made use of the equality: r (q) — Q — 1, which holds for positive q 
and for positive singularities of Holder exponents less than unity ^,18]. In turbulence, the 



most probable Holder exponent is 0.33 (corresponding to the K41 value) and for all practical 
purposes the values of Holder exponents lie between and 1 (see jlslQ)- Hence the use of 
the above equality is well justified. 

Equation El could be used to validate our previous claim, that the parameters of j^Q] give 
rise to a monoaffine field (i.e., ( q is a linear function of q). If we consider | c?x I — \d 2 \ — d — 
2- 1 / 3 , then, C, = l-log 2 {\d x \ q + \d 2 \ q ) = l-log 2 (2d«) = -glog 2 (d) = -glog 2 (2~V3) = q /3 
[QED]. Equation El could also be used to derive the classic result of Barnsley regarding 

that the graph dimension (or box- 



the fractal dimension of IFS. It is well-known 



21 



counting dimension) is related to £i as follows: D — 2 — £i- Now, using Equation El we get, 



A? 



D = 2 — Ci = 1 + l°gjv X] \dn\- F° r N = 2, we recover Equation 01 



n=l 



Intuitively, by prescribing several scaling exponents, ( q (which are known apriori from 
observational data), it is possible to solve for d n from the overdetermined system of equa- 
tions (Equation El). These solved parameters, d n , along with other easily derivable (from 
the given anchor points and d n ) parameters (a n , c n , e n and f n ) in turn can be used to 
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construct multiaffine signals. For example, solving for the values quoted by Frisch 
C 2 = 0.70, C 3 = 1, C 4 = 1-28, Cs = 1-53, Ce = 1-77, (7 = 2.01 and Cs = 2.23, along with 
£1 = 0.33 (corresponding to D — 5/3), yields the stretching factors \d n \ = 0.887,0.676. 
There are altogether eight possible sign combinations for the above stretching parameter 
magnitudes and all of them can potentially produce multiaffine fields with the aforemen- 
tioned scaling exponents. However, all of them might not be the "right" candidate from the 
LES-performance perspective. Rigorous a-priori and a-posteriori testing of these Multiaffine 
SGS models is needed to elucidate this issue (see Section |VJ). 

We repeated our previous numerical experiment with the stretching parameters d\ = 
—0.887 and d 2 = 0.676. Figure 3a shows the measured values (ensemble averaged over one 
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hundred realizations) of the scaling exponents ( q upto 12 th order. For comparison we have 
also shown the theoretical values computed directly from Equation El (dashed line). A model 
proposed by She and Leveque [23] based on a hierarchy of fluctuation structures associated 
with the vortex filaments is also shown for comparison (dotted line) . We chose this particular 
model because of its remarkable agreement with experimental data. The She and Leveque 
model predicts: Cg = f + 2 — 2(|) g / 3 . Figure 3b shows the pdfs of the increments, which are 
quite similar to what is observed in real turbulence - for large r the pdf is near Gaussian 
while for smaller r it becomes more and more peaked at the core with high tails (see also 
Figure 7b for the variation of flatness factors of the pdfs of increments with distance r). 



IV. FRACTAL CALCULUS AND SUBGRID-SCALE MODELING 



In the case of an incompressible fluid, the spatially filtered Navier-Stokes equations are: 

dltr, 



dXr, 



0. 



(7a) 



dUr 



dt 



+ u n 



dUn 

OXr. 







dXr 



-0„ 



+ r„ 



+ v\J 2 u m 



(7b) 



m, n = 1, 2, 3. 



where t is time, x n is the spatial coordinate in the n-direction, u n is the velocity component in 
the n-direction, p is the dynamic pressure, p is the density and v is the molecular viscosity 
of the fluid. The tilde ( ) denotes the filtering operation, using a filter of characteristic 



width A 



24j . These filtered equations are now amenable to numerical solution (LES) on a 



grid with mesh-size of order A, considerably larger than the smallest scale of motion (the 
Kolmogorov scale). However, the SGS stress tensor r mn in Equation [ZEj defined as 



Um.Ur 



Um.Ur. 



(8) 



is not known. It essentially represents the contribution of unresolved scales (smaller than 
A) to the total momentum transport and must be parameterized (via a SGS model) as a 
function of the resolved velocity field. Due to strong influence of the SGS parameteriza- 
tions on the dynamics of the resolved turbulence, considerable research efforts have been 



made during the past decades and several SGS models have been proposed (see [26|, |2jj for 
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reviews). The Eddy- viscosity model (13) and its variants (e.g., the Dynamic model j^], 
the Scale-Dependent Dynamic model [30]) are perhaps the most widely used SGS models. 
They parameterize the SGS stresses as being proportional to the resolved velocity gradi- 
ents. These SGS models and other standard models (e.g., Similarity, Nonlinear, Mixed 
models) postulate the form of the SGS stress tensors rather than the structure of the SGS 
fields ((31 1). Philosophically a very different approach would be to explicitly reconstruct the 
subgrid-scales from given resolved scale values (by exploiting the statistical structures of the 
unresolved turbulent fields) using a specific mathematical tool (e.g., the Fractal Interpola- 
tion Technique) and subsequently estimate the relevant SGS tensors necessary for LES. The 
Fractal model of (| ?| and our proposed Multiaffine model basically represent this new class 



26 



27 



m 



of SGS modeling, also known as the "direct modeling of SGS turbulence" 

In Section Hm we have demonstrated that FIT could be effectively used to generate syn- 
thetic turbulence fields with desirable statistical properties. In addition, Barnsley's rigorous 
fractal calculus offers the ability to analytically evaluate any statistical moment of these 
synthetically generated fields, which in turn could be used for SGS modeling. Detailed 
discussion of the fractal calculus is beyond the scope of this paper. Below, we briefly sum- 
marize the equations most relevant to the present work. Let us first consider the moment 
i 

integral: U^ m = f x m [u{x)} dx. In the present context (xo = 0,xi = 0.5 and X2 = 1), 
o 

this moment integral could be viewed as 2A filtering (A = 0.5) with the Top Hat filter 
[ i.e., Fa {x) = 1/A, if \x\ < A/2 and F&(x) = 0, otherwise ]. For instance, the 1-D com- 
ponent of the SGS stress tensor reads as: 

t = uu —uu (9a) 
= E/2,0 -#1,0^1,0 (9b) 

Barnsley (jsj|) proved that for the fractal interpolation IFS (Equation [TJ, the moment 
integral becomes: 



Lm 



m—1 I Tfl \ 2 '~ 1 l+m~ p 

E U ld £ «i +1 4C" j + E E K(l,m,p,j)U PJ 

j=0 \ j I n=1 P=° 3=0 



1 - £ <HX 

n=l 



(10a) 



s 



where, 

2 / , \ l+m-p 
^ I \a n (c n x + f n ) l ~~ P dl (a n x + e n ) m = ^ K (I, m,p,j) x 1 (10b) 

n=l \P J j=0 

After some algebraic manipulations, the SGS stress equation at node Xi becomes: 

Ti = «0Mi_l + OClUi + Oi 2 uf +1 + 

a 3 Mi_i£ti + a 4 UiU i+1 + a 5 u i+1 Ui-i (11) 

We would like to point out that the coefficients a& are sole functions of the stretching factors 
d n . In other words, if one can specify the values of d n in advance, the SGS stress (r) could 
be explicitly written in terms of the coarse-grained (resolved) velocity field (v,i) weighted 
according to weights uniquely determined by d n . In Table HI we have listed the values 
corresponding to eight stretching factor combinations, \d n \ = 0.887,0.676. It is evident 
that any two combinations (g^,^) and (o?2,^i) are simply "mirror" images of each other in 
terms of Thus, only four distinct Multiaffme SGS models (Ml, M2, M3 and M4) could 
be formed from the aforementioned eight \d n \ combinations and in each case the orderings 
could be chosen at random with equal probabilities. In this table, we have also included 



the fractal model of 0, Q| and the Similarity model of 



33| in expanded form similar to 



the Multiaffme models (see the Appendix for more information on standard SGS models). 
The Multiaffme models and the Fractal model differ slightly in terms of filtering operation, 
performed filtering at a scale A (see Equation IA.4aJ) . whereas in the case of Similarity 



model, |34j found that it is more appropriate to filter at 2 A. For the Multiaffme models, we 
also chose to employ 2A filtering scale. 

One noticable feature in Table U is that some combinations of d n result in strongly asym- 
metric weights a.}.. As an example, in the case of M4 with d\ = +0.676 and c?2 = +0.887, 
|«o| ^> local an d \ocs\ ^> |ct 4 1 . This means that the SGS stress at any node X{ would have 
more weight from the resolved velocity at node than node Xj+i. One would expect that 
such an asymmetry could have serious implication in terms of SGS model performance. 

fn the following section, we will attempt to address this issue among others by evaluating 
several SGS models via the a-priori analysis approach. 
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TABLE I: The Multiaffine, Fractal and Similarity SGS models in expanded form and their corre- 
sponding coefficients for the computation of SGS stresses according to Equation II II 



Model 




(di, 


d 2 ) 


Filter Width 


0:q 








0'2 


«3 


«4 


«5 


Multiaffine (Ml) 


(-0 


887, 


+0 


676) 


2A 





218 





204 





050 


— 


372 


— 


036 


-0 


065 




(+0 


676, 


-0 


887) 


2A 





050 





204 





218 


-0 


036 


-0 


372 


-0 


065 


Multiaffine (M2) 


(+0 


887, 


-0 


676) 


2A 





030 





248 





261 


-0 


018 


-0 


479 


-0 


043 




(-0 


676, 


+0 


887) 


2A 





261 





248 





030 


— 


479 


-0 


018 


-0 


043 


Multiaffine (M3) 


(-0 


887, 


-0 


676) 


2A 


o 


144 





220 





133 


-0 


230 


-0 


209 


-0 


057 




(-0 


676, 


-0 


887) 


2A 





133 





220 





144 


-0 


209 


-0 


230 


-0 


057 


Multiaffine (M4) 


(+0 


887, 


+0 


676) 


2A 





064 





319 





262 


-0 


121 


-0 


517 


-0 


007 




(+0 


676, 


+0 


887) 


2A 





262 





319 





064 


-0 


517 


-0 


121 


-0 


007 


Fractal 


(-0 


794, 


+0 


794) 


A 





127 





221 





026 


-0 


322 


-0 


120 


+0 


069 




(+0 


794, 


-0 


794) 


A 





026 





221 





127 


-0 


120 


-0 


322 


+0 


069 


Similarity 




NA 




2A 





188 





250 





188 


-0 


250 


-0 


250 


-0 


125 



V. EVALUATION OF SGS MODELS: A-PRIORI ANALYSIS APPROACH 



The SGS models and their underlying hypotheses can be evaluated by two approaches: 
a- priori testing and a-posteriori testing (terms coined by |35|). In a-posteriori testing, LES 
computations are actually performed with proposed SGS models and validated against refer- 
ence solutions (in terms of mean velocity, scalar and stress distributions, spectra etc.). How- 
ever, owing to the multitude of factors involved in any numerical simulation (e.g., numerical 
discretizations, time integrations, averaging and filtering), a-posteriori tests in general do 
not provide much insight about the detailed physics of the newly implemented SGS models 
26, 27 1 . A complementary and perhaps more fundamental approach 2(| would be to use 



high-resolution model (Direct Numerical Simulation or DNS), experimental or field observa- 
tional data to compute the "real" and modeled SGS tensors directly from their definitions 
and compare them subsequently. This approach, widely known as the a-priori analysis, does 
not require any actual LES modeling and is theoretically more tractable. In this work we 
focused on comparing the performance of SGS models via the a-priori analysis. We strictly 
followed the 1-D apriori analysis approach of To highlight the caveats of the 
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proposed and several existing SGS models, we performed an extensive inter-model compar- 
ison study. This exercise also helped selecting the "right" combination of stretching factors 
for the Multiaffine SGS models. 

In general, the correlation between real (r Real ) and modeled (r Model ) SGS stresses is 
considered to be a good indicator of the expected performance of a proposed SGS model. 
Another crucial indicator is the so called SGS energy dissipation rate (II): 

II = - Tl3 S l3 (12) 

15 9m ^ . . . 

t— (1-D approximation) 

2 ox 

In the inertial range, the SGS energy dissipation rate is the most influential factor affecting 
the dynamical evolution of the resolved kinetic energy |26j. On average, II is positive, 
representing a net drain of resolved kinetic energy into unresolved motion. Intermittent 
negative values of IT, known as "backscatter" , implies energy transfer from SGS to resolved 
scales. Unfortunately, a high correlation between real and modeled SGS stress (or SGS 
energy dissipation rate) is not a sufficient condition for the success of a proposed LES SGS 
model, although it is a highly desirable feature (27, 34, 36]. 

We primarily made use of an extensive atmospheric boundary layer (ABL) turbulence 
dataset (comprising of fast-response sonic anemometer data) collected by various researchers 
from the Johns Hopkins University, the University of California-Davis and the University 
of Iowa during Davis 1994, 1995, 1996, 1999 and Iowa 1998 field studies. Comprehensive 
description of these field experiments (e.g., surface cover, fetch, instrumentation, sampling 



frequency) can be found in 



We further augmented this dataset with nocturnal ABL 



turbulence data from CASES-99 (Cooperative Atmosphere- Surface Exchange Study 1999), 



3- 



For 



a cooperative field campaign conducted near Leon, Kansas during October 1999 
our analyses four levels (1.5, 5, 10 and 20m) of sonic anemometer data from the 60m tower 
and the adjacent mini-tower collected during two intensive observational periods (nights of 
October 17 th and 19 th ) were considered [the sonic anemometer at 1.5m was moved to 0.5m 
level on October 19 th ] . Briefly, the collective attributes of the field dataset explored in this 
study are as follows: (i) surface cover: bare soil, grass and beans; (ii) sampling frequency: 
18 to 60 Hz; (iii) sampling period: 20 to 30 minutes; (iv) sensor height (z): 0.5 to 20m; and 
(v) atmospheric stability (z/L, L is the local Obukhov length): (neutral) to 10 (very 
stable). 
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ABL field measurements are seldom free from mesoscale disturbances, wave activities, 
nonstationarities etc. The situation could be further aggravated by several kinds of sensor 
errors (e.g., random spikes, amplitude resolution error, drop outs, discontinuities etc.). Thus, 
stringent quality control and preprocessing of field data is of utmost importance for any rig- 
orous statistical analysis. Our quality control and preprocessing strategies are qualitatively 



similar to the suggestions of 



4l|. Specifically, we follow these steps: 



(1) Visual inspection of individual data series for detection of spikes, amplitude resolution 
error, drop outs and discontinuities. Discard suspected data series from further analyses. 

(2) Adjust for changes in wind direction by aligning sonic anemometer data using 60 
seconds local averages of the longitudinal and transverse component of velocity. 

(3) Partitioning of turbulent-mesoscale motion using discrete wavelet transform 
(Symmlet-8 wavelet) with a gap scale of 100 seconds. 

(4) Finally, to check for nonstationarities of the partitioned series, we performed the 
following step: we subdivided each series in 6 equal intervals and computed the standard 
deviation of each sub-series (oj, i — 1 : 6). Hmax(ai)/min(<Ji) > 2, the series was discarded. 

After all these quality control and preprocessing steps, we were left with only 358 "re- 
liable" series for a-priori analyses. These streamwise velocity series were filtered with a 
Top-Hat filter (A = 1,2,4 or 8m) and downsampled on the scale of the LES grid (A) to ob- 
tain the resolved velocity field Ui In a similar way, the streamwise SGS stress, r Rea \ was 
computed from its definition (Equation 9a). Filtering operations were always performed in 
time and interpreted as 1-D spatial filtering in the streamwise direction by means of Taylor's 
frozen flow hypothesis. The spatial derivatives were also computed from the time deriva- 
tives by invoking Taylor's hypothesis: = — where, < u > is the mean streamwise 
velocity. 

In Figure 4 representative realizations of the real and several modeled SGS stresses are 
presented. The modeled SGS stresses, r Model were computed from the definitions given in 
the Appendix. Along the same lines, the real and modeled SGS energy dissipation rates 
(Figure 5) were calculated according to U Real = -f r^ff and U Model = - f T MoM 
respectively. The SGS model constants like C$ of the Smagorinsky model or Cl of the 
Similarity model (see Appendix) were obtained by matching the mean real and modeled 

]. For consistency, the same procedure was followed 



SGS energy dissipation rates |36, 



37 



for the SGS-kinetic-energy-based model, Fractal model and Multiafnne models. In other 
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words, we always ensured that < -^r Kea '|f > = < -}^ T Modeim > Note that? thig 
procedure has no effect on the correlation results presented below. In actual simulations 
these model coefficients could be obtained dynamically following the approach of [29. \^\. 

From Figures 4 and 5 it is visually evident that both Similarity and Multiaffine models 
capture the variability of the SGS stress and energy dissipation rates reasonably well. On 
the other hand, the performances of the Smagorinsky and SGS-kinetic-energy-based models 
are very poor. Note that, the Smagorinsky model assumes that the trace of the SGS tensor 
is subtracted from the tensor, which is not feasible in 1-D a-priori analysis [3| • Thus, direct 
magnitude-wise comparison between real and the Smagorinsky model based SGS stress or 
dissipation energy is not possible. However, this does not prevent us from quantifying the 
performance by correlation coefficient. Moreover, the Smagorinsky model is by construction 
fully dissipative. Hence, this model is unable to reproduce the backscatter effects (see Figure 
5b), which do occur in the real SGS dissipation series (Figure 5a). 

In Table |HJ for A = lm, we show the correlation between the real and modeled SGS 
stress and energy dissipation rates. The standard deviations are given in parentheses. The 
model M3 is significantly better than any other Multiaffine model and this could only be 
attributed to its near-symmetric stencil structure (see Table 1). This resolves our previous 
dilemma regarding the selection of one Multiaffine SGS model from a class of four. From 
here on, we'll only report results for M3 and will identify it as the Multiaffine model. 

Next, in Figure 6, we plot the mean correlation between real and modeled SGS stress 
and energy dissipation rates for A = 1,2,4 and 8m. As anticipated, for all the models, 
the correlation decreases with increasing filtering scale. Also, the correlation of real vs. 
model SGS energy dissipation rates is usually higher compared to the SGS stress scenario, 
as noticed by other researchers. 

It is expected that in the ABL the scaling exponent values (( q ) would deviate from the 
values reported in [ljj due to near- wall effect. This means that the stretching factors d n 
based on the ( q values we used in this work, are possibly in error. Nevertheless, the overall 
performance of the Multiaffine model is beyond our expectations. It remains to be seen how 
the proposed SGS scheme will perform in a-posteriori analysis and such work is currently in 
progress. 



13 



TABLE II: Average correlation between observed and modeled SGS stresses and energy dissipation 
rates (A = lm). The results are based on 358 ABL turbulent velocity series measured during several 
field campaigns. The quantities in the parenthesis represent standard deviation. 

Corr(7~-R ea ' rj-Model^ Corr(II^ ea ' Yl^°^ e ^ 

Smagorinsky 0.25(0.09) 0.41(0.17) 

Similarity 0.49(0.10) 0.76(0.15) 

SGS-KE 0.23(0.08) 0.42(0.17) 

Fractal 0.33(0.05) 0.61(0.06) 

Multiaffine (Ml) 0.44(0.05) 0.71(0.05) 

Multiaffine (M2) 0.40(0.05) 0.68(0.06) 

Multiaffine (M3) 0.49(0.05) 0.77(0.05) 

Multiaffine (M4) 0.42(0.05) 0.70(0.05) 



VI. PASSIVE SCALAR 

Our scheme could be easily extended to synthetic passive-scalar (any diffusive component 
in a fluid flow that has no dynamical effect on the fluid motion itself, e.g., a pollutant in 
air, temperature in a weakly heated flow, a dye mixed in a turbulent jet or moisture mixing 
m air 

HQ) 

field generation. The statistical and dynamical characteristics (anisotropy, 
intermittency, pdfs etc.) of passive-scalars are surprisingly different from the underlying 



) ot pas 

0, Q 



turbulent velocity field |45J, M 

\. For example, it is even possible 



or the passive-scalar 



field to exhibit intermittency in a purely Gaussian velocity field Similar to the 

K41, neglecting intermittency, the Kolmogorov-Obukhov-Corrsin (KOC) hypothesis predicts 
that at high Reynolds and Peclet numbers, the g^-order passive-scalar structure function 
will behave as: (\6 (x + r) — 6 (x)\ 9 ) ~ ri in the inertial range. Experimental observations 
reveal that analogous to turbulent velocity, passive-scalars also exhibit anomalous scaling 
(departure from the KOC scaling). Observational data also suggest that passive-scalar fields 

nn 

are much more intermittent than velocity fields and result in stronger anomaly 45, 46]. 

To generate synthetic passive-scalar fields, we need to determine the stretching parameters 
d\ and d 2 from prescribed scaling exponents, ( q . Unlike the velocity scaling exponents, the 
published values (based on experimental observations) of higher-order passive-scalar scaling 
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exponents display significant scatter. Thus for our purpose, we used the predictions of a 

I — I 

newly proposed passive-scalar model 47]: C? = 2 + (|) 2 — 2 (f)^ 6 — (|) 2 (j^Y^- This 



model based on the hierarchical structure theory of 23] shows reasonable agreement with 
the observed data. Moreover, unlike other models, this model manages to predict that the 
scaling exponent ( q is a nondecreasing function of q. Theoretically, this is crucial because, 
otherwise, if ( q — > — oo as q — > +oo, the passive- 

scalar field cannot be bounded QH- 
Employing Equation (jUJ) and the scaling exponents (upto 8 th -order) predicted by the 
above model, we get the following stretching factors: \d n \ = 0.964, 0.606. We again repeated 
the numerical experiment of section IIHI and selected the stretching parameter combination: 
d\ = —0.964 and d 2 = 0.606. Like before, we compared the estimated [using Equation (4)] 
scaling exponents from one hundred realizations with the theoretical values [from Equation 
(6)] and the agreement was found to be highly satisfactory. To check whether a generated 
passive-scalar field (d\ = —0.964, d 2 = 0.606) possesses more non-Gaussian characteristics 
than its velocity counterpart (d\ = —0.887, d 2 = 0.676), we performed a simple numerical 
experiment. We generated both the velocity and passive-scalar fields from identical anchor 
points and also computed the corresponding flatness factors, K, as a function of distance r 
(see Figure 7b). Comparing Fig 7a with Fig 3b and also from Fig 7b, one could conclude 
that the passive-scalar field exhibits stronger non-Gaussian behavior than the velocity field, 
in accord with the literature. 



VII. CONCLUDING REMARKS 

In this paper, we propose a simple yet efficient scheme to generate synthetic turbulent ve- 
locity and passive-scalar fields. This method is competitive with most of the other synthetic 



turbulence emulator schemes (e.g. 



mm 



] ) in terms of capturing small-scale properties 



of turbulence and scalars (e.g., multiaffinity and non-Gaussian characteristics of the pdf of 
velocity and scalar increments). Moreover, extensive a-priori analyses of field measurements 
unveil the fact that this scheme could be effectively used as a SGS model in LES. Poten- 
tially, the proposed Multiaffine SGS model can address two of the unresolved issues in LES: 
it can systematically account for the near-wall and atmospheric stability effects on the SGS 
dynamics. Of course, this would require some kind of universal dependence of the scaling 
exponents on both wall-normal distance and stability. Quest for this kind of universality 



15 



has began only recently |48l . |4' 
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APPENDIX 

The standard Smagorinsky eddy-viscosity model is of the form: 

T Smag = _ 2 ( Cs A) 2 \S\S ij (Ala) 

where, 

~ 1 duj duj 
11 2 [ dx j + dx i ) 

is the resolved strain rate tensor and 

\S\ = (2SijSij) 

is the magnitude of the resolved strain rate tensor. Cs is the so-called Smagorinsky coeffi- 
cient. 

For 1-D surrogate SGS stress: 

~ du 

J 11 — TT- 

OX 

Further, by assuming that the smallest scales of the resolved motion are isotropic, the 
following equality holds [oil ]: 

< SijSij >= — < sfi > 

Employing this assumption for the instantaneous fields, we can write 

\S\ = {2~S ij S ij ) 1 I^V^\^\ 
16 



Hence, the Smagorinsky SGS stress equation becomes: 



r ^ = -2(C 5 A) 2 v^5|g|(g) (A.lb) 



The second model we considered is the Similarity Model 



33 



34]. 



rt 3 imL = C L {u { u 3 - u t Uj) (A. 2a) 

The overbar denotes explicit filtering with a filter of width 7 A (usually 7 = 2). Cl is the 
Similarity model coefficient. 

The 1-D surrogate SGS stress could be simply written as: 

r Siml = C L (m-uu) (A. 2b) 

Now, for 2A filtering this equation becomes: 

T Siml _ Q L y ^iz± ~*~ ^+1 ) _ + 2ttj + Mj+l ^ 2 j ^ 2 C ) 

which on further simplification leads to the expression in Table H] 

T Sml = C L [0.188uti + 0.255? + 0.188u 2 +1 
— 0.25 , Uj_i'Uj — 0.25-UjUj+i — 0.125£tj + i£tj_i] 

Next, we consider a SGS model based on the SGS kinetic energy (q 2 ) Isi|: 

r Stke = _ 2C 2 A f^f^ (A3a) 
v OX 

where, q 2 = 3(u — u) 2 

Here, CV is the SGS model coefficient. 

In the case of the PVaeta. .node, offlQ, the nnUnown sn bgri d s,e 88 W prodnced by 
synthetic fractal field around any grid point be written as: 

4 4 

r f ™ c = I [ u (x)] 2 dx - [ / u{x)dx} 2 (A.4a) 



1 1 
4 4 



' 48 °i u5 i u + 

1 + 15df - 24df + 12d? _ 2 ^ 2 ... 
1^2(1^!) (A - 4b) 
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where, 5iU = {u i+ i — {tj_i)/2, 5fu = u i+1 — 2ui + and = ±2 x / 3 . 
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FIG. 1: (a) A synthetic turbulence series of fractal dimension D = 5/3. The black dots denote 
initial interpolating points, (b) Structure functions of order 2, 4 and 6 (as labeled) computed 
from the series in Figure la. The slopes (( q ) corresponding to this particular realization are 0.62, 
1.25 and 1.89, respectively, (c) Pdfs of the normalized increments of the series in Figure la. The 
plus signs correspond to r = 2~ 14 , while the circles refer to a distance r = 2~ 6 . The solid curve 
designates the Gaussian distribution for reference. 
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FIG. 2: Wavelet power spectrum (cirlces) of the series in Figure la. The -5/3 power law is also 
shown for comparison. 
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6 q " 10 12 11 5u/<(8u r )V' 2 

FIG. 3: (a) The scaling exponent function £ 9 . The continuous, dashed and dotted lines denote the 
K41, Equation and the She-Leveque model predictions respectively. The circles with error bars 
(one standard deviation) are estimated values over one hundred realizations using d\ = —0.887 
and d-2 = 0.676. Experimental data of Anselmet et al.'s [5] is also shown for reference (star signs), 
(b) Pdfs of the normalized increments of the multiaffine series. The plus signs denote r = 2 -14 , 
while the circles refer to a distance r = 2~ 6 . The solid curve designates the Gaussian distribution 
for reference. 
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FIG. 4: A comparison of the real and modeled SGS stresses, computed from atmospheric boundary 
layer measurements, using 1-D filtering and Taylor's hypothesis. The filter width A is 2m. (a) 
Real, (b) Smagorinsky model, (c) Similarity model, (d) SGS kinetic energy based model, (e) Fractal 
model, and (f) Multiaffine model (M3). 
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FIG. 5: A comparison of the real and modeled SGS energy dissipation rates, computed from 
atmospheric boundary layer measurements, using 1-D filtering and Taylor's hypothesis. The filter 
width A is 2m. (a) Real, (b) Smagorinsky model, (c) Similarity model, (d) SGS kinetic energy 
based model, (e) Fractal model, and (f) Multiaffine model (M3). 
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FIG. 6: (a) Correlation between observed and modeled subgrid-scale stresses and (b) correlation 
between observed and modeled subgrid-scale energy dissipations as a function of filter width A. 
The results are based on 358 ABL turbulent velocity series measured during several field campaigns. 
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FIG. 7: (a) Pdfs of the normalized increments of the passive scalar Multiaffine series. The plus 
signs refer to distance r = 2 -14 , while the circles to a distance r = 2 -6 . The solid curve designates 
the Gaussian distribution for reference, (b) The flatness factors of the pdfs of the increments of 
the velocity (circles) and passive-scalar field (stars) as a function of distance r. Note that both the 
fields approach the Gaussian value of 3 only at large separation distances. Clearly the passive-scalar 
field is more non-Gaussian than the velocity field. 
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